Introduction to ftmsRanalysis

Lisa Bramer, Amanda White, Allison Thompson, Lee Ann McCue

Introduction

FT-MS Analysis

[TODO: Lisa this needs an edit] Fourier-transform mass spectrometry (FT-MS) is a type of mass spectrometery for determining the mass-to-charge ratio (m/z) of ions based on the cyclotron frequency of the ions in a fixed magnetic field. Chemical composition can be determined for a portion of the observed peaks/mass-to-charge ratios. FT-MS instrument data can be interpreted as peak intensities for each observed peak. FT-MS analysis has been used to examine a wide range of complex mixtures, including soils, plants, aquatic samples, petroleum and various beverages.

The ftmsRanalysis package was designed to help with various steps of processing FT-MS data, including:

• data formatting and manipulation
  
• reproducible analysis pipeline
  
• filtering data based on various properties
  
• calculating meta information for each peak (e.g. nominal oxidation state of Carbon)
  
• data visualization and summary
  • one sample
    
  • multiple samples
    
  • group comparisons

Example data

An example dataset has been included with the ftmsRanalysis package. This dataset is a subset of an experiment to assess differences in soil organic matter between multiple locations and crop types. The data were collected from two locations (M and W) for two crop flora (S and C). The data were analyzed with a 12T FTICR (Fourier-transform ion cyclotron resonance) mass spectrometer.

Data loading

Experimental data

Data required for the ftmsRanalysis package is comprised of three data tables:

• Expression Data - observed data for each peak (rows) and sample (columns)
  • values of each cell represent the peak intensity observed
• Sample Data - data capturing relevant experimental factors (columns) for each sample (rows) * e.g. samples and their sampling locations, treatment applied, etc.
  
• Molecular Identification Data - other characteristics or quantified values (columns) for each peak (rows)
  • e.g. molecular formulae

e_data (Expression Data)

The edata object is a data frame with one row per peak and one column per sample. It must have one column that is a unique ID (e.g. Mass).

library(ftmsRanalysis)
data("ftms12T_edata")
str(ftms12T_edata)

## 'data.frame':    24442 obs. of  21 variables:
##  $ Mass         : num  98.5 98.8 98.8 101.7 103.3 ...
##  $ EM0011_sample: num  0 0 5524739 0 0 ...
##  $ EM0013_sample: num  0 13070372 0 0 0 ...
##  $ EM0015_sample: num  0.0 0.0 2.4e+07 0.0 0.0 ...
##  $ EM0017_sample: num  0 16120890 0 0 0 ...
##  $ EM0019_sample: num  0 21228496 0 0 0 ...
##  $ EM0061_sample: num  1197974 0 30656158 0 0 ...
##  $ EM0063_sample: num  0 12305626 0 0 0 ...
##  $ EM0065_sample: num  0.0 1.1e+07 0.0 0.0 0.0 ...
##  $ EM0067_sample: num  0 0 12664590 0 0 ...
##  $ EM0069_sample: num  2535836 38329628 0 0 0 ...
##  $ EW0111_sample: num  0 0 21416774 0 0 ...
##  $ EW0113_sample: num  0 8070914 0 0 0 ...
##  $ EW0115_sample: num  3636046 0 38608164 0 0 ...
##  $ EW0117_sample: num  0 3965230 0 0 0 ...
##  $ EW0119_sample: num  0 0 2439325 0 1153547 ...
##  $ EW0161_sample: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ EW0163_sample: num  0 0 0 0 0 0 0 0 0 0 ...
##  $ EW0165_sample: num  0 0 0 16443347 0 ...
##  $ EW0167_sample: num  0 1598118 0 0 0 ...
##  $ EW0169_sample: num  0 0 0 0 0 0 0 0 0 0 ...

f_data (Sample Data)

The fdata object is a data frame with one row per sample with information about experimental conditions. It must have a column that matches the sample column names in edata.

data("ftms12T_fdata")
str(ftms12T_fdata)

## 'data.frame':    20 obs. of  4 variables:
##  $ SampleID  : Factor w/ 20 levels "EM0011_sample",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Location  : Factor w/ 2 levels "M","W": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Block     : int  1 2 3 4 5 1 2 3 4 5 ...
##  $ Crop.Flora: Factor w/ 2 levels "C","S": 2 2 2 2 2 1 1 1 1 1 ...

e_meta (Molecular Identification Data)

The emeta object is a data frame with one row per peak and columns containing other meta data. Either a column giving the molecular formula or elemental count columns (currently the elements C, H, O, S, N, and P are supported) are required. It must have an ID column corresponding to the ID column in edata. If information about isotopic peaks is available and specified, these peaks are currently filtered from the data upon peakData object creation.

data("ftms12T_emeta")
str(ftms12T_emeta)

## 'data.frame':    24442 obs. of  10 variables:
##  $ Mass       : num  98.5 98.8 98.8 101.7 103.3 ...
##  $ C          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ H          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ O          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ N          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ C13        : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ S          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ P          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Error      : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ NeutralMass: num  99.5 99.8 99.8 102.7 104.4 ...

Constructing a peakData object

peakObj <- as.peakData(ftms12T_edata, ftms12T_fdata, ftms12T_emeta, 
                       edata_cname="Mass", fdata_cname="SampleID", 
                       mass_cname="Mass", c_cname="C", h_cname="H", 
                       o_cname="O", n_cname="N", s_cname="S", 
                       p_cname="P", isotopic_cname = "C13", 
                       isotopic_notation = "1")
peakObj

## peakData object
## # Peaks: 23060
## # Samples: 20
## Meta data columns: [Mass, C, H, O, N, C13, S, P, Error, NeutralMass, MolForm]

The as.peakData function also allows for the following (optional) parameters:

• data_scale - assumed to be ‘abundance’ or peak intensity. Other options include: log2, log10, log, presence/absence (0/1)
  
• instrument_type - assumed to be 12T/15T. The option is 21T for which data is displayed differently in visualizations due to the high resolution of the data
• extraction_cname - name of column in f_data specifying extraction (e.g. water)
• isotopic_cname - name of column in e_meta which indicates if a peak is isotopic
• isotopic_notation - character string in isotopic_cname which indicates a peak is isotopic. Isotopes are currently filtered out of the data

The resulting peakData object contains three elements, named e_data, f_data, and e_meta:

names(peakObj)

## [1] "e_data" "f_data" "e_meta"

During object construction, the molecular formula is calculated from the elemental columns (and elemental columns would be created in the case that molecular formulae were provided):

tail(peakObj$e_meta)

##              Mass  C  H  O N C13 S P     Error NeutralMass     MolForm
## 24437 897.1796269  0  0  0 0   0 0 0 0.0000000    898.1869        <NA>
## 24438 897.2209292  0  0  0 0   0 0 0 0.0000000    898.2282        <NA>
## 24439 897.3973977 36 69 22 1   0 0 1 0.2345417    898.4047 C36H69NO22P
## 24440  898.812526  0  0  0 0   0 0 0 0.0000000    899.8198        <NA>
## 24441 899.0458907  0  0  0 0   0 0 0 0.0000000    900.0532        <NA>
## 24442 899.3370941  0  0  0 0   0 0 0 0.0000000    900.3444        <NA>

There is a summary method:

summary(peakObj)

## Samples: 20
## Molecules: 23060
## Percent Missing: 81.739%

… and a default plot method:

plot(peakObj)

Preprocessing

Transforming abundance values

When dealing with ’omics data quantitatively, we often log-transform to stabilize variances and reduce skew for downstream data processing. Alternatively, it’s common to treat FT-MS data as presence/absence data. We can use the edata_transform function to transform the data scale for either of these options.

peakObj <- edata_transform(peakObj, data_scale="log2")

# for presence/absence transformation:
edata_transform(peakObj, data_scale="pres")

## Warning in edata_transform(peakObj, data_scale = "pres"): Data has already
## been transformed. Data is on the log2 scale. Changing to pres scale.

## peakData object
## # Peaks: 23060
## # Samples: 20
## Meta data columns: [Mass, C, H, O, N, C13, S, P, Error, NeutralMass, MolForm]

When we plot the transformed data, the difference in scale is apparent.

plot(peakObj)

Calculating meta-data

It is frequently useful for biological analysis and interpretation to calculate values related to chemical properties of each peak, such as the nominal oxidation state of Carbon (NOSC), aromaticity, and elemental ratios. This can be done via the compound_calcs function. By default, this function calculates all available meta-data fields, specific fields can be chosen with the calc_fns parameter.

peakObj <- compound_calcs(peakObj)
peakObj

## peakData object
## # Peaks: 23060
## # Samples: 20
## Meta data columns: [Mass, C, H, O, N, C13, S, P, Error, NeutralMass, MolForm, AI, AI_Mod, DBE, DBE_O, DBE_AI, GFE, kmass, kdefect, NOSC, OtoC_ratio, HtoC_ratio, NtoC_ratio, PtoC_ratio, NtoP_ratio]

Classification of compounds based on their elemental composition is often desirable. The assign_elemental_composition function accomplishes this task.

peakObj <- assign_elemental_composition(peakObj)
table(peakObj$e_meta[,getElCompColName(peakObj)])

## 
##   CHNO  CHNOP  CHNOS CHNOSP    CHO   CHOP   CHOS  CHOSP 
##   3728   1416    742    317   3528     95    884    249

Further, each compound formula can also be assigned to biochemical compound classes (e.g. lipids, lignins, etc.) based on their chemical properities (e.g. O:C, H:C ratios), and the assign_class function performs this assignment.

peakObj <- assign_class(peakObj, boundary_set = "bs1")
table(peakObj$e_meta[, getBS1ColName(peakObj)])

## 
##                        Amino Sugar                       Carbohydrate 
##                                352                                698 
##           Carbohydrate;Amino Sugar                   Cond Hydrocarbon 
##                                 64                               1745 
##                             Lignin                 Lignin;Amino Sugar 
##                               2901                                103 
##                      Lignin;Tannin                              Lipid 
##                                 20                                709 
##                      Lipid;Protein                              Other 
##                                 21                               1840 
##                            Protein                Protein;Amino Sugar 
##                               1206                                 35 
##                     Protein;Lignin         Protein;Lignin;Amino Sugar 
##                                  7                                  2 
##                             Tannin            Tannin;Cond Hydrocarbon 
##                                984                                 57 
##                  Unsat Hydrocarbon Unsat Hydrocarbon;Cond Hydrocarbon 
##                                203                                 12

There are three sets of class boundary definitions that may be used (for the boundary_set parameter) corresponding to the following publications:

• bs1 - Kim, S., et al (2003). Analytical Chemistry.
• bs2 - Bailey, V. et al (2017). Soil Biology and Biochemistry.
• bs3 - Rivas-Ubach, A., et al (2018). Analytical chemistry.

Filtering

There are multiple types of filtering algorithms provided in ftmsRanalysis:

• Molecule filter: filter rows of e_data to exclude rows observed in too few samples
• Mass filter: filter rows based on mass, e.g. to reflect observational sensitivity of the instrument
• Formula filter: filter rows based on whether they have a molecular formula
• Emeta filter: filter rows of e_data based on a quantity/column in e_meta

For example, to filter peaks to include only masses between 200 and 900:

filter_obj <- mass_filter(peakObj)
plot(filter_obj, min_mass=200, max_mass=900)

summary(peakObj)

## Samples: 20
## Molecules: 23060
## Percent Missing: 81.739%

peakObj <- applyFilt(filter_obj, peakObj, min_mass = 200, 
                  max_mass = 900)
summary(peakObj)

## Samples: 20
## Molecules: 19327
## Percent Missing: 79.299%

Other filtering options include number of molecule observations, formula presence or absence, or emeta columns.

peakObj <- applyFilt(molecule_filter(peakObj), peakObj, min_num=2)
peakObj <- applyFilt(formula_filter(peakObj), peakObj)
peakObj <- applyFilt(emeta_filter(peakObj, "NOSC"), peakObj, min_val = 0.5)
summary(peakObj)

## Samples: 20
## Molecules: 1521
## Percent Missing: 50.352%

Visualizations of one sample

To construct plots of a single sample, first we must subset the data to contain only one sample, via the subset method.

one_sample <- subset(peakObj, samples="EM0011_sample")
summary(one_sample)

## Samples: 1
## Molecules: 1521
## Percent Missing: 57.791%

head(one_sample$e_data)

##             Mass EM0011_sample
## 3746 200.9433045            NA
## 3748 200.9863723            NA
## 3774 202.9413892            NA
## 3844 209.0091827      20.32056
## 3892 212.0200788      20.82018
## 3909  213.004142            NA

Van Krevelen plot

A Van Krevelen plot shows H:C ratio vs O:C ratio for each peak observed in a sample that has a molecular formula (thus H:C and O:C are known). By default, the points are colored according to molecular composition class determined by bs1.

vanKrevelenPlot(one_sample, title="EM0011_sample")

By default, this function colors by Van Krevelen class. However, we can also color the points according to other meta data columns in the e_meta object.

vanKrevelenPlot(one_sample, colorCName="PtoC_ratio", 
                title="Color by P:C Ratio", legendTitle = "P:C Ratio")

Kendrick plot

A Kendrick plot shows Kendrick Defect vs Kendrick mass for each observed peak.

Ions of the same family have the same Kendrick mass defect and are positioned along a horizontal line on the plot. A Kendrick plot is often used in conjunction with a Van Krevelen plot for evaluating elemental composition.

kendrickPlot(one_sample, title="Kendrick Plot for EM0011_sample")

Histogram

We can also plot the distributions of any (numeric) column of meta-data (i.e. column of e_meta).

densityPlot(one_sample, variable = "NOSC", plot_curve=TRUE, plot_hist=TRUE,
            title="NOSC Distribution for EM0011_sample")

It’s also possible to plot just the histogram or just the density curve with this function with the plot_hist and plot_curve parameters.

densityPlot(one_sample, variable = "kmass", 
            title="Kendrick Mass for EM0011_sample", plot_hist=TRUE, 
            plot_curve = FALSE, yaxis="count")

Comparison of experimental groups

The goal of this experiment was to identify differences in soil organic matter between sample locations and crop types. In order to do that we need to compare experimental treatments (groups).

The group_designation method defines treatment groups based on the variable(s) specified as main effects. Here we define groups based on the crop/flora type.

peakObj <- group_designation(peakObj, main_effects=c("Crop.Flora"))
getGroupDF(peakObj)

##         SampleID Group
## 1  EM0011_sample     S
## 2  EM0013_sample     S
## 3  EM0015_sample     S
## 4  EM0017_sample     S
## 5  EM0019_sample     S
## 6  EM0061_sample     C
## 7  EM0063_sample     C
## 8  EM0065_sample     C
## 9  EM0067_sample     C
## 10 EM0069_sample     C
## 11 EW0111_sample     C
## 12 EW0113_sample     C
## 13 EW0115_sample     C
## 14 EW0117_sample     C
## 15 EW0119_sample     C
## 16 EW0161_sample     S
## 17 EW0163_sample     S
## 18 EW0165_sample     S
## 19 EW0167_sample     S
## 20 EW0169_sample     S

The summarizeGroups function calculates group-level summaries per peak, such as the number or proportion of samples in which peak is observed. The resulting object’s e_data element contains one column per group, per summary function.

group_summary <- summarizeGroups(peakObj, summary_functions = 
                                   c("n_present", "prop_present"))
head(group_summary$e_data)

##          Mass S_n_present S_prop_present C_n_present C_prop_present
## 1 200.9433045           3            0.3           4            0.4
## 2 200.9863723           2            0.2           0            0.0
## 3 202.9413892           1            0.1           1            0.1
## 4 209.0091827          10            1.0           7            0.7
## 5 212.0200788           9            0.9           9            0.9
## 6  213.004142           2            0.2           3            0.3

We can use the densityPlot function to compare distributions of a molecular property (e.g. NOSC) between groups.

densityPlot(peakObj, samples=FALSE, groups=c("S","C"), variable="NOSC", 
            title="Comparison of NOSC Between Crop Types") 

We might also want to look at which peaks occur only in one group or another, versus those that appear in both groups. The number or proportion of samples for which a peak must be observed can be specified to determine if a peak was present for a group. Similarly, a threshold based on the number or proportion of samples can be specified to determine when a peak is absent from a group. Alternatively, a statistical test called the G-Test can be used. This is a likelihood ratio test which tests the hypothesis that the presence/absence of a peak across samples is independent of group membership.

The first step is to create peakData objects that each contain two groups to facilitate group comparisons

byGroup <- divideByGroupComparisons(peakObj, 
                                comparisons = "all")[[1]]$value

crop_unique <- summarizeGroupComparisons(byGroup, 
            summary_functions="uniqueness_gtest", 
            summary_function_params=list(
                  uniqueness_gtest=list(pres_fn="nsamps", 
                          pres_thresh=2, pvalue_thresh=0.05)))

tail(crop_unique$e_data)

##             Mass uniqueness_gtest
## 1516 653.0299655             <NA>
## 1517 657.0215477             <NA>
## 1518 658.1682292 Observed in Both
## 1519 658.1686302 Observed in Both
## 1520 700.9749846             <NA>
## 1521  719.098301             <NA>

Then we can construct a Van Krevelen plot colored by group membership.

vanKrevelenPlot(crop_unique, colorCName = "uniqueness_gtest")

The same could be done with a Kendrick plot.